#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 18 16:03:29 2022

@author: liqingsimac
"""

####################################################
####################################################
## 9.3. 有限差分空间导数示意图
'''
import numpy as np
import matplotlib.pyplot as plt

x=np.linspace(0,1,11)
y=np.linspace(0,0.2,11)
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(x,np.zeros_like(x),'bo-')
ax.set_xlim(left=-0.2,right=1.2)
ax.set_ylim(bottom=-0.05,top=0.25)
ax.plot(np.zeros_like(y),y,'bo-')
ax.plot(1+np.zeros_like(y),y,'bo-')

for k in range(9):
    yy=y[0:4]
    xx=np.zeros_like(yy)+k*0.1+0.1
    ax.plot(xx,yy,'g-')

ax.plot(x[1:-1],np.zeros_like(x[1:-1])+0.02,'ro')
ax.plot(x,np.zeros_like(x)+0.02,'r-')
ax.set_xlabel('x = space variable')
ax.set_ylabel('t = time')
'''

####################################################
####################################################
## 9.4. 一个函数求导的例子，比较差分法与谱方法
'''
import numpy as np
from scipy.fftpack import diff

def fd(u):
    #Return 2*dx finite difference x-derivative of u.
    ud=np.empty_like(u)
    ud[1:-1]=u[2: ]-u[ :-2]
    ud[0]=u[1]-u[-1]
    ud[-1]=u[0]-u[-2]
    return ud

for N in [4,8,16,32,64,128,256]:
    dx=2.0*np.pi/N
    x=np.linspace(0,2.0*np.pi,N,endpoint=False)
    u=np.exp(np.sin(x))
    du_ex=np.cos(x)*u
    du_sp=diff(u)
    du_fd=fd(u)/(2.0*dx)
    err_sp=np.max(np.abs(du_sp-du_ex))
    err_fd=np.max(np.abs(du_fd-du_ex))
    #print('N=%d,err_sp=%.4e err_fd=%.4e'%(N,err_sp,err_fd))
    print('N=%d & %.2e & %.2e \hline \\ '%(N,err_sp,err_fd))

N=32
x=np.linspace(0,2.0*np.pi,N+1)
u=np.exp(np.sin(x))

import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(x,u,'bo-')
ax.set_xticks(np.linspace(0,2.0*np.pi,4+1))
ax.set_ylim(bottom=0,top=3.0)
'''

####################################################
####################################################
## 9.5. 对空间变量呈现周期性的偏微分方程
'''
import numpy as np
from scipy.fftpack import diff
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def u_exact(t,x):
    #Exact solution.
    return np.exp(np.sin(x-2*np.pi*t))

def rhs(u,t):
    #Return rhs.
    return -2.0*np.pi*diff(u)

N=16
x=np.linspace(0,2*np.pi,N,endpoint=False)
u0=u_exact(0,x)
t_initial=0.0
t_final=64*np.pi
t=np.linspace(t_initial,t_final,51)
sol=odeint(rhs,u0,t,mxstep=500)

u_ex=u_exact(t[-1],x)
err=np.max(np.abs(sol[-1,: ]-u_ex))
print('With %d Fourier nodes the final error = %g' %(N,err))

fig=plt.figure()
ax=Axes3D(fig)
t_gr,x_gr=np.meshgrid(x,t)
ax.plot_surface(t_gr,x_gr,sol,cmap='viridis')
ax.elev,ax.azim=60,-140
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')

fig2=plt.figure()
bx=fig2.add_subplot(111)
bx.plot(x,u_ex,'b.-',label='exact solution')
bx.plot(x,sol[-1,: ]+0.1,'r.-',label='numerical solution + 0.1')
bx.set_xlabel('x')
bx.set_ylabel('u(t_final,x)')
bx.legend()
bx.set_title('exact u(tn,x) and numerical u(tn,x)')


fig3=plt.figure()
cx=fig3.add_subplot(111)
for k in range(5):
    cx.plot(x,sol[k*2,: ],'.-')
cx.set_xlabel('x')
cx.set_ylabel('u(t,x)')
cx.set_title('u(t,x) for t=t0,t1,t2,t3,t4')
'''

####################################################
####################################################









